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1. INTOOOUCTION 


The conputer s'liulation of pi cture-e lenient (pixel) data sets in remote 
sensing problems is an important method for testing proposed oiethods of 
data analysis. In such singulations the distribution of the number of pixels 
to be assigned to a rare signatu'-e O'" crop type can be modeled by a Poisson 
process. Unfortunately, Poisson pseudo-random generators for moderate mean 
rates ("*50) tend to be either slw or ir.a 'jratt. This technical memorandum 
proposes an alternative method of generation to partially overcome these 
problems . 
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2. APPROXIMATELY POISSON PSEUOO-RANOOM VARIABLES 

X 

A Poisson-distributed random variable x with probability e X /x 
has mean and variance Because of the iT^ortance of Poisson processes, 
computational procedures for generating pseudo-random Poisson-distributed 
numbers have received some attention in the literature {Snow 1968, Schaffer 
1970). The methods oescribed tiiere work ^uue satisfactorily for A small, 
but involves a sequence of calculation- of deration order A. Thus, as A 
becomes large, the time required to qener at. ■ a single Poisson variate may be 
prohibitive. The IMSL program library (IMSL 1979) which implements the 
above procedurfc replaces them with a normally distributed random variable 
with A mean anu v;;; ince tor a >5U. As wii * be seen later, this approximation 
is not particularly good. 

The criterion for a good approximation lo a particular distribution will be 
the ■»- metric 

maxlF (x)-F(x) Iwhere F(x) is the cummulative distribution 

function of the random variable to be approximated. Since satisfactory 
procedures exist t.:, generate normal pseudo-random variables, we will use 
approximations of the following form: 

i; I Normal (0,1), then = T(p+oZ) 

where T is an appropriate one-one transformation and p and o are constants 
depending on T and A. Thus, if ’i> is the normal fO, 1) cumulative distribution 
function. 


F = 4>[{T Vx)-h)/n] 

approx ' ' ' " 

The Edgeworth expansion (Abramovitz and Stegun 1970) of the function T \x) 

(where x is Poisson distributed) shows that the highest order correction 

-1 3 ■ 

term to the normality of T (x) is the skewness Y| "= (where the 

third centra! moment) times a function of x independent of T. Thus, we wish 

to find a T such that T ^(x) tias skewness as small as possible. 
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3. POWER TRANSFORMATIONS 

The fanily of functions T ^ which we will explore will be the powers 
x^, a>o, which are all defined since x Poisson > 0. This 
choice of family is suggested by the well-known "variance-stabilizing" 
transformation^ wdose var iance is constant to order X \ This has been 
useful in analyzing Pnisson-cuunt data since after transformation the 
separate counts may be treated as hr^tnosredastic. 

This corresponds to our case a = 1/2. To estimate the skewness of x^, 

X Poisson, we write 





and using the binomial theorem 

. t ...J 


then using standard results about the central moments of a Poisson variable 


= 4 . a(g-1 ) ,a-l a(a-l)(a-2) (3a-5) ,a-2 

■ 2 ^ 24 ^ +••• 


Manipulation of this expansion gives us 

, ,a . ^a-1) ,a-l ^ 
u a - A + ^ A + . . . 

X 2 

o^a - aA*^’’/^ + 3/4a(a-l)^A^"^^^ + ... 


^Ix'^ = a(3a-2)A~’^‘' 
Thus, by choosing a = 2/3, x^ 


+ ... 

is asymptotically unskewed as X becomes large. 
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4. A GENERATION ALGORITHM 


The random number generation scheme suggested by this result is as follows: 

Generate a normal (0,1) pseudo-random variate Z. Multiply it by 

a = 2/3X^^^ + 1/18 a'^'^^ and add p = - 1/9X”^^^ 

then raise the result to the 3/2-power. The Poisson variate x is taken to be 
the non-negative integer nearest to this value. 

In table 1, the «>-metric error for this approximate scheme is tabulated 
for various values of X, along with the error for the usual untransformed 
(a = 1) normal approximation. It is apparent that the 3/2-scheme is 
enormously more accurate; for X as small as 2 it is already better than 
the usual procedure at X - 50. Also, as ' increases the 3/2-scheme becomes 
better as X , whereas the old scheme becomes better as X Prichard 

(1980) has implemented this scheme on a programmable calculator with very 
little difficulty. 

The disadvantages of the scheme include the extra computation involved in 
calculating m and a, which may be negligible if a number of variates are 
to be calculated, and the overhead involved in taking a 3/2 power for each 
variate. The higher accuracy may well compensate for this. In mixed 
exact-and-asymptotic for small and large X schemes, speed may actually be 
enhanced because the asymptotic approximation can be tolerated for much 
smaller X's. 
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TABLE 1. -COMPARISON OF MAXIMUM ERRORS IN TWO APPROXIMATIONS 
TO THE POISSON DISTRIBUTION 


Max Error In 

NortPal ( Norma 1)^^^ 

X Ap prox. Approx. 


2 

.044 

.0077 

5 

.0290 

.00293 

10 

.0208 

.00144 

15 

.0170 

.00095 

20 

.0148 

.000695 

25 

.0132 

.000556 

30 

.0121 

.000457 

40 

.0105 

.000341 

50 

.000938 

.000271 

75 

.00766 

.000178 

100 

.00664 

.000132 

150 

.00495 

.000072 
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5. EaIENSIONS 


As a further application of the results of this memorandum, we may construct 
a goodness-of-f it test for catecjories that are believed to have independent 
Poisson processes as generators. Let 0^ , i = 1 , . . . , n be the observed 
counts and E- be the expected couats (^'s). Then 


E 


i = l 



0^/3 (E^) 


2 

\ 

•'n 



where the approximation is mu-wli cIoS'T '<i'" small [ . ’s than is the Pearson 

2 ^ 

X . Linear constraints, however, are no longet linear sO this test warrants 

further investigation. 
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